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ABSTRACT 

In  [1]  and  [2],  algorithms  were  introduced  for  adaptively  computing 
smooth  piecewise  polynomial  approximations  using  uniform,  least-squares  (8.2) » 
and  restricted  range  uniform  approximations.  This  paper  deals  with  the  intro- 
duction of  adaptive  curve  fitting.  Adopting  the  least-squares  algorithm, 
a package  has  been  developed  that  will  compute  smooth  piecewise  polynomial 
approximations  to  data  and/or  precise  mathematical  functions  (in  discrete 
form)  allowing  the  user  the  option  of  using  best  t-j  or  best  ^2  approximations. 
(The  code  for  this  newly  developed  adaptive  curve  fitting  package 

supercedes  the  code  listed  in  [1]  for  the  old  least-squares  package.)  The 
adaptive  curve  fitting  algorithm  used  in  this  package  is  described  in  detail 
in  section  II.  Section  III  involves  a discussion  of  the  -approximation 
problem  as  a linear  programming  problem,  as  well  as  a description  of  the  j^i 
algorithm  used  in  the  computing  of  the  t-j  approximations.  In  section  IV, 
the  FORTRAN  program  that  has  been  developed  for  this  t]-t2  adaptive  curve 
fitting  algorithm  is  discussed,  and  numerical  results  are  presented  in  an 
effort  to  illustrate  how  the  version  of  this  algorithm  may  be  used  most 
effectively. 


I.  Introduction 

Let  X be  a finite  set  of  real  points  and  let  f be  a function  defined  on  X, 

or  let  data  be  given  in  tabular  form.  In  the  case  of  data  given  in  tabular 

form,  say  {(t^. , we  shall  let  X = and  define  f on  X by  f(t^) 

• y^»  1 ® 1.  2,  M,  so  that  the  functional  notation  may  be  used  in  what 

follows.  Let  a = min{x:  x € X}  and  b = max{x:  x €X}  and  for  any  function  g 

defined  on  X denote  max{g(x):  x € X}  by  Hgllj^.  Finally,  let  N,  SMTH,  and  TOL 

be  parameters  supplied  by  the  user  where  N and  SMTH  are  integer  values  with 

N ^ SMTH  ^ -1  and  TOL  is  a positive  number.  In  this  setting,  the  adaptive 

curve  fitting  algorithm  presented  here  will  calculate  a piecewise  polynomial 

1/ 

approximation  p to  f and  a set  of  points  (knots)  C X with  a = Xq  < x^  < 

...  < X|^  = b such  that: 

(1)  p restricted  to  [x^_p  x^.]  is  a polynomial  p^-  e = {q:  q is  a 
real  algebraic  polynomial  of  degree  ^N-l), 

(2)  p has  SMTH  continuous  derivatives,  and 

(3)  ||f  - pIIx  iTOL. 

II.  The  Adaptive  Curve  Fitting  Algorithm 

We  shall  begin  with  a somewhat  brief  description  of  the  algorithm,  in  an 
attempt  to  illustrate  the  basic  logical  flow.  This  shall  be  followed  by  an 
in-depth  discussion  of  the  algorithm's  various  components.  Again,  it  should 
be  mentioned  that  the  input  to  the  algorithm  consists  of  the  function  to  be 
approximated  (in  the  discrete  form  as  described  above),  and  the  values  N,  SMTH, 
and  TOL. 

The  algorithm  begins  by  finding  the  largest  point  x^  in  X such  that: 

(i)  [fli  x^]  nx  contains  at  least  max(2,  N+1)  points,  and 
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(11)  the  best  or  Hg)  approximation  p^  € .j  to  f on  [a,  x-j]  H X meets 
the  presdribed  tolerance  (i.e.,  11^  - P^ll^g  x-jlnx 

If  x-|  ■ b,  then  since  p^  is  a piecewise  polynomial  satisfying  (l)-(3),  the 
algorithm  is  successfully  terminated.  If  x-j  < b,  then  the  right  endpoint  of 
the  first  subinterval  (i.e.,  the  knot  x^)  is  generally  determined  by  "backing 
off"  from  x^  to  a point  in  (a,  x-jln  X selected  in  such  a manner  (to  be  explained 

later)  so  as  to  add  to  the  stability  of  the  algorithm. 

Having  thus  determined  the  first  subinterval  [a,  x-j]  and  corresponding 
best  (t^  or  ^2)  approximation  p^  € to  f on  [a,  x-j]  H X,  the  algorithm  goes 

on  to  determine  the  remainder  of  the  piecewise  polynomial  approximation  p. 

However,  for  the  remainder  of  the  approximation,  the  problem  of  meeting  the 
smoothness  constraint  SMTH  at  the  interior  knots  x-j,  Xg,  ....  X|^__^  must  be  con- 
sidered. Thus,  the  next  step  is  to  find  the  largest  point  X2  in  X such  that: 

(1)  [x^,  X2]n  X contains  at  least  max(2,  N-SMTH)  points,  and 

(11)  the  best  (t^  or  ^2)  approximation  P2  e -j  to  f on  [x^ , Xg]  A X 

subject  to  the  constraint  that  (x^)  = pj’^^(xi),  j = 0,  1,  .... 
SMTH,  satisfies  Hf  - P2ll|-^  ^ j < TOL. 

Again,  if  X2  = b,  p^  and  P2  constitute  a piecewise  polynomial  p satisfying 
(l)-(3)  and  the  algorithm  is  successfully  terminated.  If  X2  < b,  then  the  knot 
X2  Is  determined  in  the  same  manner  as  x-j,  in  order  to  establish  the  second  sub- 
interval [Xp  X2]. 

The  algorithm  continues  in  the  same  manner  to  its  completion  by  finding 
successive  subintervals  [x2,  X3],  [x^,  x^],  ...,  [X|^_i , b],  and  corresponding 
polynomial  approximations  P3,  p^,  ...,  to  f such  that  the  TOL  and 

SMTH  constraints  are  met. 

At  this  point,  a more  detailed  description  of  the  various  components  of 
the  algorithm  is  in  order,  beginning  with  how  the  x^.  are  determined.  Suppose 


M 
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some  knot  has  been  established.  (This  may  be  the  initial  knot  Xq  = 
some  interior  knot.)  The  process  of  finding  x^  is  initiated  by  computing  the 
best  or  ^2)  approximation  (subject  to  the  smoothness  constraint  SMTH  if 
x^_1  is  an  interior  knot)  to  f on  b]  A X.  If  this  approximation  satis- 

fies the  prescribed  tolerance  TOL,  then  x^  * b.  Should  this  approximation  not 
satisfy  TOL,  the  algorithm  sets  & = b.  In  what  follows,  ^ shall  denote  the 
current  smallest  point  tn  X such  that  the  appropriate  best  approximation  to  f on 
[x^_1,  &]  n X does  not  satisfy  TOL.  Next,  the  appropriate  best  approximation  to  f 
on  [Xj_p  t]  n X is  computed,  where  t = infix  € X:  [x^_i,  x]  n X contains  at 
least  max(2,  N-SMTH)  points).  If  this  approximation  fails  to  satisfy  TOL,  the 
algorithm  cannot  meet  the  desired  accuracy  and  fails.  Otherwise,  the  algorithm 
sets  a = t.  In  what  follows,  a shall  denote  the  current  largest  point  in  X 
such  that  the  appropriate  best  approximation  to  f on  [x^^-j,  a]  O X satisfies  TOL. 

An  iterative  procedure  then  ensues  to  find  the  largest  possible  x^.  The 
algorithm  sets  t = infix  € X:  (fe  - a)/2  £ x < b).  If  this  set  is  empty,  the 
algorithm  sets  t = suplx  € X:  a £ x £ (&  - a)/2}.  In  effect,  a "halfway  point" 
between  the  current  values  of  a and  is  being  sought.  If  t = a,  then  this 
procedure  has  converged  and  X.  = a.  Otherwise,  the  appropriate  best  approximation 
to  f on  [x^_p  t]  n X is  computed.  If  this  approximation  satisfies  TOL,  we 
have  a new  value  for  a,  namely  a = t.  If  this  approximation  does  not  satisfy 
TOL,  then  we  have  a new  value  for  fe,  namely  & = t.  This  process  is  continued, 
essentially  drawing  the  values  of  the  largest  "good"  candidate  for  x^.  (the 
current  a)  and  the  smallest  "bad"  candidate  for  x^  (the  current  ^)  closer 
together.  When  - a is  less  than  or  equal  to  some  user  definable  prescribed 
tolerance  (ETA),  or  when  a and  6 are  adjacent  elements  of  X (which  ever  occurs 
first),  the  process  is  completed  and  a is  accepted  as  a good  approximation  to  x 
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Remaining  in  the  setting  just  outlined,  i.e.,  having  just  established  an 
5(^,  the  method  of  "backing  off"  mentioned  earlier  used  to  determine  the  knot 
location  (and  thereby  establish  the  subinterval  x^.])  will  now  be 

explained.  (It  should  be  noted  here  that  if  SMTH  = -1  or  0,  the  following 
procedure  is  bypassed,  and  the  knot  location  x^  is  simply  the  point  x^.) 

First,  the  error  function  f(x)  - p^(x)  is  examined  for  those  points 
5it  ^2*  ^^i-1’  ^i^  n X at  which  relative  extrema  occur.  That  is, 

is  such  that  |f(5^)  - P^-(C^)|  1 K(x)  - p^(x)|  for  x = max{t  € X:  t < 5^} 
and  X = min{t  c X:  t > (Note  that  p^  is  the  appropriate  best  approxi- 

mation to  f on  [x^_i,  x^]n  X computed  in  the  previous  process.)  In  a process 
to  be  subsequently  described,  one  of  the  Cy  is  chosen  for  the  knot  location  x^. 
The  motivation  for  choosing  x^  in  this  manner  is  as  follows:  If  f were  differ- 
entiable and  X = [a,  b],  then  choosing  the  knot  x^  to  be  an  interior  extreme 
point  of  f - p^  on  [x^-_i,  x^]  would  ensure  that  f'(x^)  * p](x^-).  That  is,  the 
slope  of  the  function  f and  its  polynomial  approximation  on  Cx^_i»  x^]i  P^-. 
would  match  at  x^.  This  becomes  advantageous  when  joining  the  next  polynomial 
piece  of  the  approximation  to  p^.  at  x^-  when  the  approximation  is  required  to 
be  at  least  continuously  differentiable.  When  the  next  polynomial  piece,  P^-+] » 
is  smoothly  joined  to  p^  at  x^,  it  follows  that  the  function  f and  both  of  the 
polynomial  pieces  p^.  and  p^^^  shall  have  the  same  slope  about  the  knot  x^.  If 
we  simply  joined  the  polynomial  pieces  p^  and  p.^^  at  x^-,  this  is  generally  not 
the  case  and  severe  oscillatory  problems  tend  to  set  in.  This  process  of 
"backing  off"  from  x^  to  a smaller  x^  e ^^v^v=l  co'^^'"‘'t>utes  significantly  to 
the  stability  of  the  algorithm. 

In  determining  which  5y  should  be  chosen  for  the  knot  location  x^,  steps 
must  be  taken  to  alleviate  the  fact  that  f (in  the  discrete  form  input  to  the 
algorithm)  is  in  fact  not  differentiable.  The  values  3'^® 


calculated  where  T is  the  derivative  of  the  quadratic  interpolation  of  f <5— 

centered  at  Then,  is  chosen  to  be  the  largest  5^  such  that 

|f^(Cy)  - P|(5^)|  is  less  than  some  user  definable  prescribed  tolerance  (EPS). 

If  there  does  not  exist  such  a then  x^  is  chosen  to  be  the  largest  5 at 
which  |f^(5^)  - pj(5^)|  is  a minimum.  (It  should  be  noted  here  that  in  the 
Implementation  of  the  algorithm,  it  is  generally  not  the  case  that  all  of  the 
relative  extreme  points  of  f - p^  on  (x^_p  x^]  H X are  considered,  but  rather 
only  the  largest  N-SMTH-1  of  them. 

This  "backing  off"  procedure  has  proven  to  be  invaluable  for  dampening 
oscillations  for  those  approximation  in  which  SMTH  > 0.  However,  if  no  conti- 
nuity is  required  of  the  approximation  (SMTH  = -1)  or  if  the  approximation  is 
simply  required  to  be  continuous  (SMTH  = 0),  this  proceoure  serves  no  purpose 
and  is  overridden.  (Again  we  note  that  in  these  instances  x^  = x^,  i =1,  2,  .~.) 

Another  aspect  of  the  algorithm  that  merits  attention  is  the  process  of 
determining  the  last  subinterval  of  the  approximation.  Suppose  that  the  sub- 
intervals [a,  x^],  [xp  X2],  ....  Cx^_2.  and  corresponding  polynomial 

approximations  p-j , P2,  ...»  P^-_i  if  f have  been  determined.  If  [x^_i,  b]  n X 
contains  at  least  max(2,  N-SMTH)  points,  the  algorithm  goes  on  to  determine  x^ 
precisely  as  outlined  earlier.  However,  If  [x^_i.  b]  H X contains  fewer  than 
max(2,  N-SMTH)  points,  a process  is  initiated  that,  if  successful,  will  deter- 
mine the  last  subinterval  of  the  approximation.  In  this  instance,  the  knot 
x^,^  is  replaced  with  some  x^_i  in  (x^_2.  x^_^)  0 X chosen  such  that  [x^_].  b]nx 
will  contain  at  least  max(2,  N-SHTH)  points.  Specifically,  the  algorithm 
chooses  x^._i  to  be  that  point  in  X closest  to,  (b  - x^._2)/2  such  that: 

(0  f^i-l»  ^ X contains  at  least  max(2,  N-SMTH)  points,  and 

(ii)  the  best  (t-i  or  ^2)  approximation  p^-  € to  f on  [x^,-].  b]  n X 


subject  to  the  constraint  that  p|'^^(x^.l)  = pj^l(Xj_i).  j=0,l,  ...,  SMTH, 


t 

I It  Is  most  often  the  case  that  such  a point  Is  readily  available,  in  which 

[ case  the  algorithm  ts  successfully  terminated.  (Here,  the  last  two  subintervals 

I and  corresponding  polynomial  approximations  to  f are,  respectively,  [x^_2» 

( 

Pi-1’  ^i’^  ® satisfactory  x^_^  cannot  be  found,  the  algorithm 

is  terminated  and  an  appropriate  error  message  is  generated. 

I The  bulk  of  the  discussion  thus  far  has  been  centered  around  the  manner 

I in  which  the  knots  a = Xq,  x-j,  ...,  Xj^  = b of  the  piecewise  polynomial  approxi- 

mation p to  f are  determined.  This  is  due  to  the  fact  that  the  algorithm's 

adaptive  nature,  which  allows  it  to  calculate  the  total  number  and  location  of 

I 

knots  needed  (as  opposed  to  most  data  fitting  spline  techniques  that  require 
that  the  total  number  and  location  of  the  knots  be  specified  in  advance),  is 

I one  of  the  main  features  of  the  algorithm.  We  shall  now  direct  attention  to 

the  actual  computation  of  the  polynomial  pieces  p^ , P2,  ...»  P|^  that  comprise 
i the  approximation  p to  f on  X,  this  having  thus  far  been  taken  somewhat  for 

granted. 

Suppose  that  some  knot  x^._i  has  been  established,  and  the  procedure  that 
determines  x^  is  currently  in  progress.  Suppose  further  that  some  point 
t € [x^_i,  b]  X is  being  considered  as  a candidate  for  x^.  Thus,  the  next 
step  is  to  determine  the  appropriate  approximation  to  f on  t]  H X 

(which  we  shall  denote  by  p^).  Let  [x^_p  t]  H X = {t^  <_  t2  1 . • • £ t^}-  If  < 

t^  • a,  or  if  SMTH  = -1,  there  are  no  continuity  constraints  involved  and  the 
following  over-determined  system  is  constructed: 


’10  0 
1 (t2”t^) 


0 

(t2-ti) 

^v^i) 


N-1 


N-1 


f'll 

^2 

fit^) 

J 

-S- 

J 


1 (vt,) 
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The  solution  of  the  above  system  (in  the  or  sense)  yields  the  coefficients 

{c^}".l  of  a polynomial  in  the  form  p^(x)  = ^ c^(x  - t^)''"  . If  t^  is  an 

v*l 

Interior  knot  (in  which  case  a previous  polynomial  p^_i  has  been  established) 
and  SMTH  =*  k > -1,  the  smoothness  constraint  at  t-j  must  be  considered  and  the 
following  over-determined  system  is  constructed: 


(tg-t^) 


k+1 


(t2-ti ) 


N-n 


k+1 


N-1 


'k+2 


L ‘^N 


I r Pi-1^^1^  11 


The  solution  of  the  above  system  (in  the  t-j  or  sense)  yields  the  coeffi- 
cients ^ polynomial  of  the  form 


P.(x) 


k 

I 

j*0 


M) 


(tj 


'i-lV'-l 

~JT 


-(x  - t,)‘ 


N 

I 

v«k+2 


c (x  - t, 

V 1 


kV-1 


If  the  error  of  approximation  of  f by  p^  on  [x^_^,  t]  0 X is  greater  than  TOL, 
a smaller  candidate  for  x^.  (and  corresponding  polynomial  approximation)  is  sought. 
If  p^  satisfies  TOL,  the  next  step  depends  on  t.  If  the  process  used  for 
finding  x^  has  not  converged  yet,  a larger  candidate  for  x^  (and  corresponding 
polynomial  approximation)  is  sought.  If  the  process  has  converged,  t is 
accepted  as  x^  and  the  polynomial  approximation  p^  is  accepted  as  p^. 

III.  The  t^-Approximation  Problem 

The  general  (linear)  approximation  problem  can  be  stated  as  follows: 
Suppose  we  are  given  a real-valued  function  f (to  be  approximated)  defined  on  a 
finite  set  of  real  points  X = {x^,  X2.  ...,  x^^j}.  Let  (where  n < m)  be 

a set  of  real-valued  functions  defined  on  X.  Then,  for  any  set  of  real  numbers 
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C ■ {Cii  C2.  ....  c^}  we  can  define  a linear  approximating  function  of  the  form 
n 

A(C,  x)  = c.(>.(x).  The  z. -approximation  problem  is  to  find  a best  approxi- 

jal  J J I 

mation  A(C*,  x)  (i.e.,  a set  C*  = {c^,  c^,  ....  c*})  that  minimizes 

m 

(1)  I |f(x^)  - A(C.  x,.)( 

1*1 

It  Is  well  known  that  at  least  one  best  approximation  always  exists. 

I 

The  connection  between  the  z-j -approximation  problem  and  linear  programming 
is  well  established,  and  has  become  the  basis  for  many  of  the  more  effective 
algorithms.  In  the  adaptive  curve  fitting  package  presented  here,  the 
algorithm  used  for  determining  the  z^  approximations  was  developed  by  Barrodale 
and  Roberts  [3].  It  is  a modified  version  of  the  simplex  method  of  linear 
programming  and  empirical  evidence  indicates  that  it  is  computationally  superior 
to  any  other  known  algorithm  for  solving  the  z^ -approximation  problem  (see 
Barrodale  and  Roberts  [5]  for  comparative  results).  Unlike  many  z^  algorithms 
that  require  that  the  set  of  approximating  functions  be  linearly  inde- 

pendent  on  X,  or  satisfy  the  Haar  condition  on  X,  this  algorithm  can  be  used 
with  any  set  of  approximating  functions. 

Before  going  into  a description  of  this  algorithm,  we  note  here  that 
familiarity  with  the  standard  form  of  the  simplex  method  is  being  assumed  on 
the  part  of  the  reader  (see  the  Appendix,  or  [7J).  Also,  it  should  be  pointed 
out  that  the  description  presented  here  is  somewhat  superficial  and  that  the 
papers  by  Barrodale  and  Roberts  that  are  relevant  (see  [3]  for  references)  are 
recotrmended  for  further  details. 

For  the  general  z^ -approximation  problem  outlined  earlier,  let 
♦4  4 “ '|ii(x,-).  fi  = f(x,-),  and  define  non-negative  variables  u.,  v. , a,-  and  b. 

:■  ‘ 

m — - ..-rr-srrz::::; iji 
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j ■ 1,  2,  n.  Then  a best  approximation  to  f on  X corresponds  to  an 
optimal  solution  to  the  linear  programri’ng  problem: 

m 


(2) 


Minimize  z = I (u^  + v^) 


(3) 


subject  to 


* “1  - ’1 


i - 1.  2. 


■ 


(4) 


and  aj,  bj,  u^ , v.  ^ 0. 


The  formulation  of  the  above  correspondence  appears  in  Barrodale  and  Roberts 

[6]. 

The  main  modifications  to  the  standard  form  of  the  simplex  method  as 
presented  in  [3]  result  in  an  algorithm  that  exploits  the  above  linear 
programming  version  of  the  t-j-approximation  problem  in  such  a manner  that 
reduces  the  number  of  iterations  required  to  arrive  at  an  optimal  solution 
(i.e.,  a best  approximation).  Basically,  the  algorithm  allows  for  the 
passage  through  several  neighboring  simplex  vertices  (basic  feasible  solutions, 
extreme-point  solutions)  in  a single  iteration,  thus  reducing  the  number  of 
time-consuming  simplex  transformations  needed  to  arrive  at  an  optimal  solution 

To  begin  our  description  of  the  algorithm,  we  note  that  an  initial  basic 
feasible  solution  to  the  linear  programming  problem  (2)  is  immediately  avail- 
able. That  is,  if  we  denote  the  columns  of  the  simplex  tableau  corresponding 
to  (2)-(4)  by  Aj,  Bj,  j = 1,  2,  . . . , n,  and  , i = 1 , 2,  ....  m,  an 

initial  basis  is  provided  by  , Ug,  ....  (provided  that  each  f^  is  non- 
negative). Below  is  the  resulting  initial  simplex  tableau: 


L 


Here,  the  column  to  the  left  of  the  tableau  is  used  to  indicate  the  basis, 
which  initially  consists  of  the  vectors  U^.  The  upper  portion  of  the  tableau 
corresponds  to  the  constraints  (3),  and  the  bottom  (objective)  row  is  a repre- 
sentation of  the  objective  function  (2).  That  is.  if  we  consider  (2)  in  form 


i (-U.  - V.)  + z * 0 
i*l  ’ ^ 


then  it  follows  that 


m n 


^ i’  ^''i'  i^i^"  iv^+z  = 0 

i«l  j*l  ^ ^ i=l  ^ i=l  ^ i=i  ^ 


mm  m m m m m m 


which  is  equivalent  to 
m m 


^ ^^|/2.i^®2  •••  ^J/n.iK  - ^.^/l.i^h  - U *2.i^‘’2  - 

(5) 

‘ ^''l  ' ^''2  " ■ ^''m  ^ 

Thus  we  see  that  the  objective  row  of  the  initial  simplex  tableau  is  a repre- 
sentation of  the  objective  function  as  shown  in  (5).  If  any  of  the  f^.  are 
negative,  the  sign  of  the  corresponding  row  is  changed,  and  is  replaced  in 
the  basis  by  V^.  Note  also  that  Aj  = -Bj,  j *=  1,  2,  ....  n,  U.  = -V^., 
i ■ 1.  2,  ....  m,  the  sum  of  the  entries  in  the  objective  row  corresponding  to 
each  pair  of  Aj  and  Bj  is  0,  and  the  sum  of  the  entries  in  the  objective  row 
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corresponding  to  each  pair  of  and  is  -2.  This  is  utilized  by  the  FORTRAN 
program  of  this  algorithm  which  economizes  on  storage  by  using  a condensed 
form  of  the  simplex  tableau  (see  [4]). 

The  algorithm  consists  of  two  distinct  stages.  Stage  1,  which  begins 
upon  the  establishment  of  the  initial  simplex  tableau,  restricts  the  choice 
of  the  pivotal  columns  during  the  first  n iterations  to  the  vectors  A.  and 
The  vector  to  enter  the  basis  is  chosen  to  be  that  which  has  the  largest  non- 
negative entry  in  the  objective  row.  (We  note  two  things  here:  Choosing  a 
vector  to  enter  the  basis  is  equivalent  to  selecting  a variable  to  enter  the 
set  of  basic  variables.  Also,  with  regard  to  the  choice  of  vectors  to  enter 
the  basis  (i.e.,  variables  to  enter  the  set  of  basic  variables)  in  Stage  1, 
choosing  from  among  those  with  non-negative  entries  in  the  objective  row 
clearly  serves  to  decrease  the  value  of  the  objective  function,  given  the 
formulation  of  the  objective  row.)  The  vector  to  leave  the  basis  (or  equi- 
valently, the  variable  to  enter  the  set  of  non-basic  variables)  in  Stage  1 is 
chosen  from  among  the  basic  vectors  (and  V^)  by  selecting  that  vector  which 
causes  the  maximum  reduction  in  the  objective  function. 


At  the  end  of  Stage  1,  the  resulting  simplex  tableau  represents  an  approxi- 
mation that  interpolates  at  least  K of  the  data  points,  where  K is  the  number 
of  vectors  A.  or  B.  in  the  basis.  (This  results  from  the  fact  that  K of  the 

w J 

vectors  (or  V^-)  have  been  removed  from  the  basis.)  If  the  approximation 
determined  by  the  simplex  tableau  of  the  end  of  Stage  1 interpolates  more  than 
K points,  this  is  an  indication  that  the  tableau  is  degenerate.  (This  does  not 
cause  any  problems  in  practice,  however.) 


During  Stage  2,  the  non-basic  and  vectors  are  interchanged  with  the 
basic  U-  and  Vj.  The  basic  and  B.  vectors  are  not  allowed  to  leave  the 

W sj 

basis  during  this  stage.  At  each  iteration,  the  vector  chosen  to  enter  the 
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basis  Is  that  (among  the  non-basic  and  V^)  with  the  most  positive  entry  in 
the  objective  row,  and  the  vector  to  leave  the  basis  is  again  chosen  to  be  that 
(among  the  basic  and  V^)  which  causes  the  maximum  reduction  in  the  objective 
function.  The  algorithm  terminates  when  all  of  the  entries  in  the  objective 
row  are  non-positive.  (Although  this  normally  occurs  in  Stage  2,  should  it 
occur  In  Stage  1,  the  second  stage  becomes  unneccesary. ) 

Each  vector  that  enters  the  basis  during  Stage  2 determines  a data  point 
to  be  dropped  from  the  interpolating  set  while  the  corresponding  vector  leaving 
the  basis  determines  a new  point  of  interpolation.  Thus  assuming  nondegeneracy, 
the  number  of  data  points  interpolated  by  the  approximation  resulting  from 
Stage  1 is  preserved.  The  final  simplex  tableau  at  the  end  of  Stage  2 may 
contain  basic  vectors  A.  or  B.  that  have  negative  values  associated  with  them, 

J w 

thus  resulting  in  an  infeasible  solution.  This  solution  is  made  feasible  (and 
f hence  optimal)  by  interchanging  such  basic  vectors  A^  (or  B^)  with  the  corres- 

ponding  non-basic  vectors  B.  (or  A.).  Given  the  nature  of  the  A^  and  B^  vectors, 
this  Is  a simple  matter  of  changing  the  sign  of  the  appropriate  row. 

It  Is  Important  to  note  here  the  main  modification  to  the  standard  simplex 
method  that  this  algorithm  introduces.  For  any  given  iteration  (in  both  Stage 
1 and  Stage  2),  upon  having  established  the  vector  to  enter  the  basis  (i.e., 
the  pivot  column),  the  vector  chosen  to  leave  the  basis  (i.e.,  the  pivot  row) 

Is  not  in  general  that  vector  which  the  standard  simplex  method  dictates  should 
be  removed,  but  rather  is  that  vector  (from  among  the  basic  and  V^-  vectors) 
which  causes  the  maximum  reduction  in  the  objective  function.  Geometrically, 
the  simplex  transformation  which  then  follows  is  equivalent  to  a passage 
through  several  neighboring  extreme  points  of  the  feasible  region.  The 
standard  procedure  for  determining  the  vector  to  leave  the  basis  (see  the 
Appendix  or  [7])  is  modified  as  follows: 


J 
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I 

I First,  a tentative  pivot  row  Is  determined  using  the  standard  procedure. 

I If  subtracting  twice  the  value  of  the  resulting  pivot  from  the  entry  in  the 

I objective  row  of  the  pivot  column  yields  a nonpositive  result,  then  a normal 

I simplex  transformation  is  performed  on  the  tableau.  Otherwise,  twice  the  pivot 

I row  Is  subtracted  from  the  objective  row,  the  pivot  row  is  multiplied  by  -1, 

and  the  vector  (or  V^)  in  the  basis  corresponding  to  the  pivot  row  is 

! replaced  in  the  basis  by  the  corresponding  (or  U^).  This  serves  to  decrease 

i 

I the  value  of  the  objective  function  and  change  the  sign  of  the  pivot.  The 

f standard  procedure  for  finding  a pivot  row  is  applied  again,  resulting  in  a new 

( tentative  pivot.  The  previous  procedure  is  applied  until  a pivot  is  chosen 

(which  cannot  be  rejected  (i.e.,  a pivot  such  that  twice  its  value  subtracted 
from  the  entry  in  the  objective  row  of  the  pivot  column  is  nonpositive).  A 
I normal  simplex  transformation  is  then  performed  with  this  pivot. 

I 

Remark 

f In  practice,  errors  due  to  roundoff  can  result  in  a non-basic  vector  being 

chosen  as  a pivot  column,  even  though  there  are  no  positive  entries  in  the 
vector  (and  hence  no  candidate  for  a pivot).  This  is  usually  due  to  a loss  of 
significance  that  can  occur  during  simplex  transformations  on  a tableau  con- 
taining entries  differing  greatly  in  magnitude.  In  the  FORTRAN  version  of 
this  algorithm  [4],  a small  positive  tolerance  (TOLER)  is  defined  below  which 
the  magnitude  of  any  quantity  is  considered  to  be  zero.  The  algorithm  is 
terminated  prematurely  if  a pivot  column  is  encountered  that  contains  no  candi- 
dates for  a pivot  that  exceeds  the  value  TOLER.  This  occurence  is  rare,  and 
generally  indicates  that  the  current  solution  is  close  to  the  actual  optimal 
solution.  Thus,  the  program  outputs  the  current  solution  as  a reasonable 
approximation  to  the  actual  optimal  solution  in  these  instances. 
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IV.  Numerical  Results 

This  i'yi-2  ®^3ptive  curve  fitting  algorithm  has  been  implemented  as  a i 

FORTRAN  program  (standard  ASCII  code)  and  has  been  tested  on  Colorado  State 
University's  CDC  6400  and  CDC  Cyber  172.  A fully  documented  listing  of  this 
package  and  complete  discussion  of  the  program's  major  components  shall  appear 
In  a future  paper.  Following  is  a general  discussion  of  input/output  a*  well 
as  some  numerical  results. 

Aside  from  the  function  (or  data)  f to  be  approximated  (in  discrete  form), 
those  values  that  must  be  supplied  by  the  user  include:  lOPT,  an  integer 
value  (1  or  2)  specifying  the  choice  of  approximation  (best  or  best  12)',  i 

N,  the  number  of  coefficients  of  each  polynomial  piece;  and  TOL,  the  desired 
error  tolerance.  In  circumstances  where  relatively  few  data  points  are 
available,  the  program  fills  in  the  gaps  between  the  data  points  by  discretizing 

a piecewise  linear  interpolation  of  the  original  data.  Should  this  be  necessary,  ! 

i 

the  number  of  additional  "data  points"  to  be  inserted  in  this  manner  (NPTS)  i 

1 

should  also  be  supplied  by  the  user.  ' 

i 

The  output  consists  of  the  coefficients  of  the  polynomial  pieces  (where  | 

the  polynomials  are  in  standard  fom),  the  knot  locations  of  the  piecewise  j 

I 

polynomial  approximation,  and  the  error  of  approximation  for  each  polynomial  j 

j 

piece.  Also,  all  of  the  appropriate  information  about  the  piecewise  polynomial  i 

i 

approximation  is  stored  so  as  to  allow  the  user  to  evaluate  the  approximation 
at  any  point  (in  the  interval  of  approximation)  using  an  available  function 
subprogram. 

Two  other  values  used  in  the  program  that  are  not  input  but  are  user 
definable  are  the  values  ETA  and  EPS.  Recall  (section  II)  that  in  the  iterative 
procedure  used  to  determine  the  location  of  a particular  x^. , the  values  a (the 
current  largest  "good"  candidate  for  x^)  and  & (the  current  smallest  "bad" 

I 
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candidate  for  are  drawn  together,  and  when  B-a  is  less  than  or  equal  to 
the  user  prescribed  tolerance  ETA,  or  when  B and  a become  adjacent  elements  of 
the  data  set  (whichever  occurs  first),  a is  accepted  as  the  location  for  x^. 

In  the  examples  that  follow,  ETA  was  set  to  be  less  than  or  equal  to  the  mesh 
size  of  the  data,  thus  ensuring  that  the  iterative  procedure  will  continue 
until  a and  B are  in  fact  adjacent  and  the  resulting  x^  is  as  large  as  possible. 
In  the  "backing  off"  procedure  used  to  determine  the  actual  knot  locations  x^ 
(see  section  II),  x^.  is  chosen  to  be  the  largest  at  which  |f^(5y)-P|(C^)  1 
Is  less  than  the  user-prescribed  tolerance  EPS  (where  the  are  the  points 
In  at  which  the  relative  extrema  of  the  error  function  f(x)  - p-j(x) 

occur).  Empirical  evidence  indicates  that  the  value  of  .05  for  EPS  yields 
desirable  results. 

We  now  turn  to  some  examples.  Using  the  adaptive  curve  fitting  program 
with  best  approximations  (lOPT  = 1),  the  function  v'x  on  [0,  2]  was  approxi- 

mated on  201  equally  spaced  points  with  N = 6,  SMTH  = 2,  and  TOL  = .01.  Since 
this  function  is  difficult  to  approximate  by  polynomials  near  x = 0,  the 
algorithm's  ability  to  automatically  decrease  the  length  of  the  subintervals 
near  x = 0 and  then  recover  by  lengthening  them  away  from  the  origin  is  illus- 
trated. The  results  appear  below: 

Knot  locations  0.0  .06  .18  .41  .84  1.49  2.0 

Subintervals  7 pts.  13  pts.  24  pts.  44  pts.  66  pts.  52  pts. 


(It  should  be  noted  here  that  the  purpose  of  the  above  example  is  simply  to 
Illustrate  the  adaptive  nature  of  the  program  and  that  in  general,  the 
version  of  this  adaptive  curve  fitting  algorithm  is  not  particularly  suited  for 
approximating  precise  mathematical  functions.) 
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We  now  turn  our  attention  to  some  examples  that  better  illustrate  how  the 
option  of  this  adaptive  curve  fitting  program  is  most  effectively  utilized. 

As  one  might  suspect,  given  the  nature  of  approximation,  the  piecewise 

polynomial  approximations  have  shown  to  be  very  effective  for  approximating 

on  data  sets  that  contain  points  that  are  very  inaccurate  with  respect  to  the 

overall  accuracy  of  the  data.  For  the  example  in  figure  1,  the  function 

|s1n(x)|  on  [0,  Z-n]  was  discretized  (101  points)  and  "noise"  was  generated 

using  a random  number  generator.  (The  number  of  "bad"  points  and  a bound  on 

the  deviation  was  set,  but  the  location  and  magnitude  [within  the  bound]  of  the  | 

noise  was  random.)  The  pertinent  input  values  appear  at  the  top  of  the  graph  I 

and  the  execution  time  (in  CPU  seconds)  appears  below.  (In  all  the  examples  | 

that  follow,  the  data  points  are  denoted  by  X's.)  Note  that  the  noise  has 

virtually  no  ill  effect  on  the  piecewise  polynomial  Z-j  approximation.  That 

is,  the  curve  essentially  ignores  the  "bad"  points  and  passes  through  the 

"good"  points.  ! 

In  figures  2 and  3,  the  function  e’'  on  [0,  2]  was  discretized  (101  points), 

1 

noise  was  introduced  in  the  manner  just  described,  and  the  resulting  data  set 
was  approximated  on  by  the  adaptive  curve  fitting  program,  exercising  both  the 
and  Zg  (least-squares)  option.  In  so  much  as  least  squares  approximations 
have  the  tendency  to  dampen  the  effect  of  randomly  distributed  noise,  we  might 
expect  (and  in  fact  do  achieve)  desirable  results  in  both  instances.  However, 
note  that  in  figure  3 the  approximation  was  visibly  affected  by  the  noise  in 


k i 


the  interval  (.4,  .8)  and  by  the  "bad"  point  near  1.0,  whereas  in  figure  2,  the 
approximation  again  is  unaffected  by  the  noise. 

We  have  also  had  some  success  in  fitting  smaller  data  sets  with  piecewise 
polynomial  approximations.  The  experimental  data  in  figures  4 and  5 involves 
the  bitumen  yield  and  gas  and  oil  yield  (as  a function  of  time)  from  oil  shale 
heated  to  a constant  temperature.  Since  relatively  few  data  points  were 
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available  In  each  data  sample  (13-20  points),  we  filled  In  the  gaps  between  the 
data  points  with  200  equally  spaced  points  by  the  manner  previously  described 
(input  NPTS  = 200).  In  the  example  in  figure  4,  TOL  was  set  to  a rather 
liberal  5.0  so  as  to  allow  the  approximation  the  freedom  of  passing  through  the 
scattering  of  data  points,  rather  than  forcing  an  interpclatory  type  fit  that 
would  have  resulted  if  the  tolerance  were  set  smaller.  For  the  example  in 
figure  5,  TOL  was  set  to  2.0,  which  still  gives  the  algorithm  a great  deal  of 
freedom  but  results  in  a fit  that  follows  the  data  more  closely. 

Remark 

It  should  be  noted  here  that  the  version  of  this  adaptive  curve 
fitting  algorithm  should  be  used  with  some  discretion.  For  data  containing 
noise  (particularly  large  data  sets),  it  has  proved  to  be  a very  effective 
approximating  scheme.  In  general,  however,  it  should  not  be  used  for  approxi- 
mating precise  mathematical  functions,  or  for  approximating  "good"  data  sets 
for  which  an  interpolatory  type  fit  would  be  desirable. 

Also,  a future  paper  is  anticipated  in  which  all  of  the  adaptive  curve 
fitting  programs  that  have  been  developed  thus  far  (restricted  range,  , 

and  shall  be  compared  in  an  effort  to  offer  some  insight  into  how  they 
may  be  used  most  efficiently. 
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A.  The  General  Linear  Programminq  Problem 

The  general  linear  programing  problem  may  be  stated  as  follows: 
Optimize  (i.e.,  maximize  or  minimize)  the  linear  form 
(A.l)  Z - C,X,  ♦ CjXj  + ...  * c„x„ 

subject  to  the  constraints 
(A. 2)  Xj  = 0,  j = 1,  2,  , n 

and  S'j'jX'j  + 3-j2^2  * — ^ *ln^n  — ^ ^1 

*21*1  * ®22*2  * «2n’‘n  i (i)  *>2 


Ami*!  * *•••■*  ^mn^n  i W »n. 

where  the  a^j,  b^  and  Cj  are  given  constants. 

The  linear  form  (A.l)  to  be  optimized  is  called  the  objective  function. 

Here  we  assume  that  at  least  one  coefficient  a-z  is  nonzero  in  each  row  and 

• J 

that  every  variable  appears  in  some  nontrivial  inequality  with  a nonzero 

coefficient.  We  make  no  special  assumptions  about  the  b^,  they  may  be 

positive,  negative,  or  zero.  It  should  also  be  pointed  out  here  that  the 

minimum  of  an  objective  function  occurs  at  the  same  set  of  values  as  the 

«•» 

maximum  of  the  negative  of  that  objective  function.  Thus,  regardless  of 
;|  the  constraints,  the  problem  of  minimizing  z = k^x^  + k2X2  + ...  + k^^x^  is 

I equivalent  to  maximizing  -z  = -k^x^  - k2X2  - ...  Also,  any  inequality 

I ^il’'!  •••  -^"in'ni^ 

) is  equivalent  to 


I -a^^x.,  -a^2’'2-  •"  -^'n"ni-^- 

;i;  A feasible  solution  to  the  linear  prooramminq  problem  is  a vector 

!i>  ’ 

ill  , X ■ (xi , X9,  . . . , x„)  that  satisfies  conditions  (A. 2)  and  (A. 3).  The  set  of 

r?[  , \ c n 
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all  feasible  solutions  to  the  linear  programming  problem  is  a convex  region  in 
n-diraensional  Euclidean  space  and  is  referred  to  as  the  feasible  region, 
or  region  of  feasibility.  The  feasible  region  can  either  be  void  (In  which 
case  no  solution  to  the  problem  exists),  a convex  polyhedron,  or  a convex 
region  which  is  unbounded  in  some  direction. 


Example  1 
Maximize 
Subject  to 
and 


Example  2 
Maximize 
Subject  to 
and 


z = 3x  + 5y 
X > 0,  y 
X + 3y  ^ 9 
X + y > 5 
2x  + y > 6 


By  introducing  non-negative  slack  variables 
express  the  problem  (A.1)-(A.3)  in  the  equivalent  form: 
Optimize  (i.e.,  maximize  or  minimize)  the  linear  form 
(A.4)  z = c,x,  * * ...  * c„x„  » 0x„^,  * ...  * Ox 

subject  to  the  constraints 
(A. 5)  x^  = 1 ^ 


. .,  x„.„,  we  can 
n+m 


n+m 


0,  j = 1,  2, 


• • « a 


n+m 
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and 

(A.6) 


“iri 

*21 * *22^2 


^ *2n"n  ^ V2 


*inl^l  * *m2^2 


mn  n 


+ (-)  » b„ 

n+ffl  m 


A basic  solution  to  (A.6)  is  a solution  obtained  by  setting  n variables 
equal  to  zero  and  solving  for  the  remaining  m variables  (provided  the  solution 
Is  unique).  The  m variables  are  called  basic  variables;  the  n that  are  pre- 
specified as  zero  are  called  non-basic  variables.  A basic  feasible  solution 
Is  a basic  solution  which  also  satisfies  (A. 5);  i.e.,  all  the  basic  variables 
are  non-negative.  An  optimal  feasible  solution  is  a feasible  solution  which 
also  maximizes  or  minimizes  (A. 4). 

With  the  assumption  that  a particular  linear  programming  problem  possesses 
an  optimal  solution,  the  following  holds: 

1.  There  is  an  extreme  point  (corner  point)  of  the  feasible  region  at 
which  the  objective  function  takes  on  its  maximum  or  minimum. 

2.  Every  basic  feasible  solution  corresponds  to  an  extreme  point  of  the 
feasible  region. 

From  the  above,  we  see  that  it  is  only  necessary  to  examine  extreme-point 
solutions  (i.e.,  basic  feasible  solutions).  Since  there  are  at  most  (JJ|)  of 
them,  we  have  an  upper  bound  to  the  number  of  possible  solutions  to  the 
problem  (A.4)-(A.6).  However,  for  large  n and  m,  evaluating  all  the  possible 
solutions  to  find  one  that  maximizes  or  minimizes  the  objective  function  is  an 
unreasonable  way  to  proceed.  The  simplex  method,  devised  by  G.  B.  Dantzig, 
is  a computational  scheme  that  selects,  in  an  orderly  fashion,  a small  subset 
of  the  possible  solutions  that  converges  to  an  optimal  solution. 

This  algorithm  is  a method  for  moving  from  one  extreme  point  to  another 
In  such  a way  that  the  objective  function  is  always  improved,  or  at  the  least. 


j 


I 


1 
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the  same.  In  a finite  number  of  steps,  (usually  between  m and  2m)  an  optimal 
feasible  solution  (If  one  exists)  Is  found.  The  simplex  method  also  Indicates 
the  existence  of  alternate  optimal  solutions,  empty  feasible  regions,  and 

i 

unbounded  solutions.  It  Is  an  extremely  effective  tool  for  solving  any  linear  j 

programming  problem.  It  should  be  noted,  however,  that  the  proof  of  Its 
effectiveness  Is  a result  of  empirical  evidence,  rather  than  underlying 
theory. 

B.  The  Simplex  Method 

It  should  be  noted  at  the  outset  that  the  manners  In  which  the  standard  I 

simplex  method  can  be  described  are  as  diverse  as  the  forms  that  linear 
programming  problems  can  take  on.  In  so  much  as  most  elementary  treatises 
on  the  simplex  method  deal  with  the  "general"  linear  programming  problem:  | 


Maximize 

z * C^X^  + C2X2  + . 

. . + C_X_ 

n n 

(B.l) 

Subject  to 

Xj  i 0,  j = 1,  2,  . 

. . , n 

(B.2) 

and 

^11^1  ^ ^12^2  ^ ••• 

* «1n*r,  i '>1 

^21^1  ^ ^22^2  • • • 

• • 

• • 

* “2n*n  i •’2 
« • 

• • 

• • 

(B.3) 

®ml^l  * V2  * 

+ a__x„  < b_ 
mn  n — m 

or.  Introducing  non 

-negative  slack  variables,  the  equality  form  of  the  general 

problem: 

Maximize 

Z = CiX-j  + C2X2  + . 

••  C-X„  + 0x„.,  + ...  + 0x„ 
n n n+1  m 

(B.4) 

Subject  to 

Xj  ^ 0,  j = 1 , 2,  . 

. . , m + n 

(B.5) 

and 

®n^l  ^ ^ 

* »ln*r.  * Vl  ■ "l 

*21^1  ^ *22^2  ^ 

• • 

* *2r’‘n  * ^2  " '’2 

• • 

(B.6) 

• • 

«ml^l  ^ ®m2’'2  " ••• 

• • 

+ a X + X = b . 

mn  n n+m  m 
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[■ 

i 
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Me  shall  proceed  along  these  lines.  (Note  that  since  any  linear  programming 
problem  can  be  expressed  In  the  above  form.  It  Is  In  some  sense  a general 
form.) 

We  begin  by  Illustrating  the  simplex  method  with  a specific  linear 
programming  problem,  this  to  be  followed  by  a general  description  applicable 
to  any  problem  of  the  form  (B.4)-(B.6).  Consider  the  problem  In  example  1 of 
Appendix  A.  By  Introducing  slack  variables,  r,  s,  and  t,  this  problem  can  be 
expressed  as: 


Naximize 

z * 3x  + 2y  + 

Or  + Os 

+ ot 

(B.7) 

Subject  to 

X ^ 0,  y ^ 0, 

r ^ 0,  ! 

5 2^  0,  t 0 

(B.8) 

and 

X + y + r 

s 

7 

X + 2y  + 

s « 

10 

(8.9) 

2x  + y 

+ t » 

12  . 

• 

Note  that  a basic  feasible  solution  is  immediately  available.  That  is,  if 
we  let  X * y * 0,  then  It  follows  that  r * 7,  s = 10  and  t = 12.  This  can  be 
represented  in  tableau  form  as 


(B.10) 


Here,  the  column  to  the  left  of  the  tableau  Is  used  to  Indicate  the  set  of 
basic  variables  (or  basis)  which  In  this  instance  consists  of  the  variables 
r,  s and  t.  The  upper  portion  of  the  tableau  represents  the  constraints  (B.9) 
and  the  lower  portion  of  the  tableau,  known  as  the  objective  row,  is  a repre- 
sentation of  the  objective  function  (B.7)  expressed  In  the  form:  -3x  - 2y  + 2 = 0 
(where  we  need  not  keep  a column  for  the  z). 


X y r s t 


1 

1 

1 

0 

0 

7 

1 

2 

0 

1 

0 

10 

© 1 

0 

0 

1 

12 

-3  -2 

0 

0 

0 

0 

29 


Hence*  at  this  stage  we  have  our  first  basic  feasible  solution:  x » 0, 
y ■ 0,  r » 7,  s * 10,  and  t = 12.  The  entry  In  the  lower  right-hand  corner 
of  (B.IO)  Indicates  that  z = 0.  Geometrically,  this  basic  feasible  solution 
corresponds  to  the  extreme  point  of  the  feasible  region  (see  example  1, 

Appendix  A)  occuring  at  the  origin.  The  simplex  method  consists  of  employing 
Gauss-Jordan  elimination  to  proceed  to  another  basic  feasible  solution  (I.e., 
to  another  corner  point  of  the  feasible  region)  In  such  a way  that  the  value 
of  the  objective  function  (z)  shall  be  Increased  (or  at  worst  remain  the  same.) 

We  see  that  given  the  way  In  which  the  objective  function  Is  expressed  In 
the  objective  row,  z can  be  increased  by  Increasing  any  variable  with  a 
negative  coefficient  In  the  objective  row.  Those  variables  having  negative 
coefficients  In  the  objective  row  are  x and  y.  Since  the  greatest  Increase 
In  z will  clearly  result  from  Increasing  that  variable  with  the  most  negative 
coefficient,  we  choose  to  Increase  the  variable  x.  The  corresponding  column 
of  the  tableau,  in  this  Instance  the  first.  Is  called  the  pivot  column. 

We  now  consider  to  what  extent  can  the  variable  x be  Increased  (we  wish 
to  Increase  x as  much  as  possible  without  violating  any  of  the  constraints). 

If  we  consider  the  equations  In  (B.9),  or  equivalently  the  tableau  given  by 
(B.IO): 

r • 7 - X 
s » 10  - X - 
t « 12  - 2x 

We  see  that  (keeping  in  mind  that  y is  presently  equal  to  0)  If  we  wish  to 


incre  'e  x,  we  can  not  allow  it  be  greater  than  6.  This  is  the  smallest  of 
the  ratios  * T’  “ X’  ®3  * Increase  x beyond  6,  then  the  last 

of  the  above  equations  dictates  that  t will  become  negative,  which  is  not 
permitted  by  (B.8).  We  find  the  e-ratios  by  dividing  the  entries  in  the 
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right-hand  column  of  the  tableau  by  the  corresponding  positive  entries  in 

the  pivot  column.  The  pivot  row  is  that  which  yields  the  smallest  non-  | 

negative  ratio;  in  this  instance,  the  third  row  of  tableau  (B.IO). 

The  entry  in  the  tableau  located  in  both  the  pivot  column  and  the  pivot  | 

row  (circled  in  (B.IO))  is  called  the  pivot.  The  first  "simplex  transformation"  | 

1 

consists  of  applying  Gauss-Jordan  elimination  to  the  tableau  in  (B.IO)  using  | 

the  given  pivot.  The  resulting  tableau  appears  below:  | 

I 

I 

i 

i 

i 

(B.n)  I 

I 

i 


X 

y 

r 

s 

t 

r 

0 

© 

1 

0 

-1/2 

1 

s 

0 

3/2 

0 

1 

-1/2 

4 

X 

1 

1/2 

0 

0 

1/2 

6 

0 

-1/2 

0 

0 

3/2 

18 

Let  us  reflect  upon  what  has  transpired.  The  variable  x,  previously  equal 
to  0,  has  been  increased  to  6.  Thus  x has  become  a basic  variable,  or 
equivalently,  has  "entered  the  basis".  The  variable  t,  previously  equal  to 
12,  has  been  decreased  to  0.  Thus  t has  become  a non-basic  variable,  or 
equivalently,  has  "left  the  basis".  The  new  tableau  represents  the  basic 
feasible  solution:  x=6,  y=0,  r=l,s*4  and  t = 0.  The  entry  in  the 
lower  right-hand  corner  of  (B.ll)  indicates  that  z has  been  increased  to  18. 

Geometrically,  the  above  simplex  transformation  corresponds  to  a transfer  to 

an  adjacent  extreme  point  of  the  feasible  region,  and  the  basic  feasible  i 

solution  yielded  by  the  tableau  (B.ll)  corresponds  to  that  extreme  point  of 

the  feasible  region  located  at  the  point  (6,  0). 

Since  one  of  the  entries  in  the  object  row  of  (B.ll)  is  negative,  it  ' 

1 

J 

follows  that  z can  be  increased  further.  The  pivot  column  for  the  next 

simplex  transformation  shall  be  the  second  column.  (This  indicates  that  y 

1 4 8 

will  enter  the  basis.)  Upon  examining  the  e-ratios  (e^  = * 2,  02  *-3-  * j, 

7 7 


I 


« 
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and  $2  12)  we  see  that  the  pivot  row  is  the  first,  corresponding  to  the 


smallest  e-ratio.  (This  indicates  that  r will  leave  the  basis.)  Performing 
a second  simplex  transformation  (using  the  pivot  circled  in  (B.ll))  yields: 


X 

r 

s 

t 

y 

0 

1 

2 

0 

-1 

2 

s 

0 

0 

-3 

1 

1 

2 

X 

1 

0 

-1 

0 

1 

5 

0 

0 

1 

0 

1 

19 

(B.12) 


Since  there  are  no  negative  entries  in  the  objective  row,  we  have  arrived 
at  the  optimal  solution.  The  final  tableau  represents  the  optimal  feasible 
solution:  x = 5,  y = 2,  r = 0,  s = 2,  and  t = 0.  The  maximal  value  of  the 
objective  function  is  19.  This  solution  corresponds  to  the  extreme  point  of 
feasible  region  located  at  the  point  (5,  2). 

At  this  point,  we  summarize  the  simplex  method  as  illustrated  in  the 
previous  example.  If  we  consider  the  general  linear  programming  problem  as 


posed 

in  (B. 

4)-(B.6),  we  see  that  an 

initial  simplex 

tableau  is 

given  by: 

• 

’‘j  • 

• • 

’'n 

’'n+l 

"n+2 

•••  Vm 

• 

Vl 

®11 

*12  ••• 

*ij 

• • • 

*ln 

1 

0 

...  0 

^+2 

• 

• 

«21 

• 

d 22  * • • 

• 

*2j 

• « • 

*2n 

• 

0 

• 

1 

• 

0 

• • 

*n+i 

• 

• 

®il 

• 

• 

• 

®i2  • • ' 

• 

• 

• • • 

*in 

• 

• 

0 

* • • 

• 

1 ...  6 
• 

• 

• 

• 

(B.13) 

*n+m 

«ml 

• 

®m2  • • • 

*mj 

• • • 

*mn 

0 

0 

...  *1 

-"1 

*C2  ••• 

'h 

• • • 

0 

0 

0 

0 

We  assume  that  Xj  ^ 0 for  j = 1 

. 2, 

• • • * 

"n+m 

and  that  b.  0 for  i = 

w 

1 f2fa  • • flTIa 

Then  a 

basic 

feasible  solution 

is  immediately 

available,  namely 

’ "j 

= 0 for 

j • 1. 

2,  .. 

. , n and 

for  i 

= 1. 

2,  . 

• • y ni  • 

We  then  proceed  as 

follows; 
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1.  Determine  whether  the  current  basic  feasible  solution  is  optimal,  i.e., 
whether  all  the  entries  in  the  objective  row  are  greater  than  or  equal  to 
zero. 

2.  If  there  is  at  least  one  negative  entry  in  the  objective  row,  choose  the 
pivot  column  to  be  that  with  the  most  negative  entry  in  the  objective  row. 
(Assume  this  to  be  the  j^^  column.)  In  the  case  of  a tie,  choose 


arbitrarily  from  among  those  tied. 

3.  Determine  the  e-ratios  (e.  = - — , i * 1,  2,  ...,  m)  for  the  positive 

J *ij 

entries  in  the  pivot  column. 


4.  Select  the  pivot  row  to  be  that  with  the  smallest  non-negative  e-ratio. 
(Assume  this  to  be  the  i^^  row.)  In  the  case  of  a tie,  choose  arbitrarily 
from  among  those  tied.  (Theoretically,  a tie  in  this  step  may  cause 
problems.  See  the  remarks  that  follow.) 

5.  Carry  out  the  simplex  transformation  by  employing  Gauss-Jordan  elimination 
with  the  pivot  a^^. 

6.  Return  to  step  1. 


Each  iteration  produces  a new  basic  feasible  solution.  The  procedure 
terminates  when  no  suitable  pivot  column  or  pivot  row  can  be  found.  If  there 
is  no  pivot  column  (i.e.,  no  negative  entry  in  the  objective  row),  then  the 
current  solution  is  optimal.  If  there  is  a suitable  pivot  column,  but  all  of 
its  entries  are  either  zero  or  negative,  this  is  an  indication  that  the 
solution  is  unbounded. 

Remark 


Note  that  performing  a simplex  transformation  when  the  choice  of  pivot 
row  is  made  arbitrarily  as  a result  of  a tie  among  potential  pivot  rows  (see 
(4)  above)  results  in  a basic  variable  (that  corresponding  to  the  row  not 


chosen  to  be  the  pivot  row)  turning  to  zero.  The  resulting  basic  feasible 

» 

solution  (which  may  or  may  not  be  optimal)  is  said  to  be  degenerate.  As 
degeneracy  is  almost  never  a computational  problem  (except  in  rare  instances 
when  endless  cycl inq  around  a loop  of  non-optimal  degenerate  basic  feasible 
solutions  occurs),  degenerate  solutions  are  treated  as  any  other  while 
computing  an  optimal  solution  via  the  simplex  method.  For  a more  complete 
discussion  of  degeneracy  and  cycling,  see  [7]  or  any  other  comparable  text 
on  linear  programming. 
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fitting  scheme  employed  is  discussed  in  detail  as  is  the  algorithm  used  for 
determining  the  Lu  aoproximations.  Finally,  numerical  results  are  oresented  in 
an  attempt  to  illustrate  the  effective  use  of  this  newl  developed  Li  curve 
fitting  package.  ^ ^ ' 


J 


